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Abstract: The origin of highest energy cosmic rays (UHECR) is yet unknown. In order 
to understand their propagation we determine the probabihty that an ultrahigh energy 
(above 5 ■ 10^^ eV) proton created at a distance r with energy E arrives at earth above 
a threshold Ec- The clustering of ultrahigh energy cosmic rays suggests that they might 
be emitted by compact sources. A statistical analysis on the source density based on 
the multiplicities is presented. The ultrahigh energy cosmic ray spectrum is consistent 
with the decay of GUT scale particles. Alternatively, we consider the possibility that a 
large fraction of the ultrahigh energy cosmic rays are decay products of Z bosons which 
were produced in the scattering of ultrahigh energy cosmic neutrinos on cosmological relic 
neutrinos. Based on this scenario we determine the required mass of the heaviest relic 
neutrino. The required ultrahigh energy neutrino flux should be detected in the near 
future by experiments such as AMANDA, RICE or the Pierre Auger Observatory. 



1. Introduction 

The interaction of protons with the microwave background predicts a drop in the cosmic ray 
flux above the GZK cutoff ~5-10^^ eV. Observations do not confirm this GZK scenario. 
Instead, about 20 events above 10^'' eV were observed by a number of experiments. Since 
above the GZK energy the attenuation length of particles is a few tens of megaparsecs if an 
ultrahigh energy cosmic ray (UHECR) is observed on earth it was most probably produced 
in our vicinity. 

Any UHECR scenario needs the knowledge of the propagation of the particles through 
the cosmic microwave background (CMB). Section 2 reports on a study about propa- 
gation and determines the probability -P(r, E, Ec) that protons created at distance r with 
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energy E reach earth above a threshold Ec- Using this P one can give the observed spec- 
trum by one numerical integration for any injection spectrum. 

It is an interesting phenomenon that the UHECR events are clustered. Usually it is 
assumed that at these high energies the galactic and extragalactic magnetic fields do not 
affect the orbit of the cosmic rays, thus they should point back to their origin within a few 
degrees. In contrast to the low energy cosmic rays one can use UHECRs for point-source 
search astronomy. Recently, a statistical analysis j2| based on the multiplicities of the 
clustered events estimated the source density. In Section 3 we we review our extention [1] 
on the above analysis. Our analytical approach gives the event clustering probabilities for 
any space, intensity and energy distribution of the sources by using a single additional 
propagation function P{r, E; Ec). 

In Section 4 we review our work ^ on the scenario that the UHECRs are coming from 
decaying superheavy particles (SP) and on the determination of their masses mx by an 
analysis of the observed UHECR spectrum. Interestingly enough mx is compatible with 
the GUT scale. 

The existence of a background gas of free relic neutrinos is predicted by cosmology. 
Ultrahigh energy neutrinos (UHEz^) scatter on relic neutrinos (Ri^) producing Z bosons, 
which can decay hadronically (Z-burst) Z^. In Section 5 we summarize our comparison of 
the predicted proton spectrum with observations and review the determination of the mass 
of the heaviest Kv via a maximum likelihood analysis. 

The details of the presented results and a more complete reference list can be found 
in Uliini. 

2. Propagation of UHECR protons 

Using pion production as the dominant effect of energy loss for protons at energies > 10^^ eV, 
ref. [7j calculated P(r, E, Ec) for three threshold energies. We extended jlj the results of [7j. 
The inelasticity of Bethe-Heitler pair production is small (~ 10~^), thus we used a contin- 
uous energy loss approximation for this process. The inelasticity of pion-photoproduction 
is larger (~ 0.2 — 0.5) in the energy range of interest, thus there are only a few tens of such 
interactions during the propagation. Due to the Poisson statistics and the spread of the 
inelasticity, we will see a spread in the energy spectrum even if the injected spectrum is 
mono-energetic. 

In our simulation protons are propagated in small steps (10 kpc), and after each step 
the energy losses due to pair production, pion production and the adiabatic expansion 
are calculated. During the simulation we keep track of the current energy of the proton 
and its total displacement. We used the following type of parametrization P(r, E, Ec) = 
exp [—a • (r/1 Mpc)'']. Fig. ^ shows a{E/Ec) and h{E/Ec) for a range of three orders of 
magnitude and for five different Ec- Just using the functions of a{E/Ec) and b{E/Ec), 
thus a parametrization of P{r, E, Ec), one can obtain the observed energy spectrum for 
any injection spectrum without additional Monte-Carlo simulation. 
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Since P{r, Ep; E) is of universal usage, we have decided to make the latest numerical 
data for the probability distribution {—)dP{r, Ep] E) /dE available for the public via the 



The propagation function can be similarly determined for photons [S] , though the nec- 
essary CPU power is approximately 300 times more than for protons. Therefore, we used 
the stochastic method to test a few cases. Usually the continuous energy loss approxima- 
tion was used. In this approximation, the energy (and number) of the detected photons 
is a unique function of the initial energy and distance, and statistical fluctuations are 
neglected. The processes that are taken into account are pair production on the diffuse ex- 
tragalactic photon background, double pair production and inverse Compton scattering of 
the produced pairs. The energy attenuation length of the photons due to these processes is 
strongly influenced by the poorely known universal radio and infrared backgrounds. These 
uncertainties influences some of our results. 

A full simulation of the photon propagation function with all the statistical fluctuations 
will be the subject of a later work. 

3. Density of sources 

The arrival directions of the UHECRs measured by experiments show some peculiar cluster- 
ing: some events are grouped within ~ 3°, the typical angular resolution of an experiment. 
Above 4 • 10^^ eV 92 cosmic ray events were detected, including 7 doublets and 2 triplets. 
Above 10^^ eV, one doublet out of 14 events were found |H]. The chance probability of 
such a clustering from uniform distribution is rather small [31^1 ■ 

The clustered features of the events initiated an interesting statistical analysis assuming 
compact UHECR sources [2j. The authors found a large number, ~ 400 for the number of 
sources within the GZK sphere. We generalize their analysis The most probable value 
for the source density is really large; however, the statistical significance of this result is 
rather weak. 

The number of UHECRs emitted by a source of luminosity A during a period T follows 
the Poisson distribution. However, not all emitted UHECRs will be detected. They might 
loose their energy during propagation or can simply go to the wrong direction. For UHE- 
CRs the energy loss is dominated by the pion production in interaction with the cosmic 
microwave background radiation. It can be taken into account with the help of the above 
mentioned probability function P{r,E,Ec). 

The features of the Poisson distribution enforce us to take into account the fact that 
the sky is not isotropically observed. 

The probability of detecting k events from a source at distance r with energy E can be 
obtained by simply including the factor P[r, E , Ec)Ar] / {A7:r'^) in the Poisson distribution: 



o 
o 
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where we introduced j = XT A/ (Air) and Ar] / (Airr'^) , which is the probabihty that an 
emitted UHECR points to a detector of area A. The factor r/ represents the visibihty of 
the source, which was determined by spherical astronomy. We denote the space, energy 
and luminosity distributions of the sources by p(x), c{E) and h{j), respectively. The 
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probability of detecting k events above the threshold Ec from a single source randomly 
positioned within a sphere of radius R is 

Pk= dV p(x) / dE c{E) / dj h{j) X 
Jsr j Ec Jo 

exp\—P(r,E,Ec)'ni/r'^] r , . . oik , . 
^ ^ [P(r, E, E,)7jj/r^] . (3.2) 

Denote the total number of sources within the sphere of sufficiently large radius (e.g. 
several times the GZK radius) by N and the number of sources that gave k detected events 
by Nk- Clearly, = Yl'o' -^i the total number of detected events is A'^e = ScT^^*- 
The probability that for N sources the number of different detected multiplets are is: 

oo ^ 

P{N,{N,}) = ml[—P^^. (3.3) 
fc=0 ^' 

For a given set of unclustered and clustered events (A'^i and N2,N^,-.-) inverting the 
P{N, {Nk}) distribution function gives the most probable value for the number 

Note, that Pk and then P{N, {Nk}) are easily determined by a well behaved four- 
dimensional numerical integration for any c{E), h{j) and p{r) distribution functions. In 
order to illustrate the uncertainties and sensitivities of the results we used a few different 
choices for these distribution functions. 

For c{E) we studied three possibilities. The most straightforward choice is the ex- 
trapolation of the 'conventional high energy component' oc E~'^. Another possibility is to 
use a stronger fall-off of the spectrum at energies just below the GZK cutoff, e.g. oc E^^^. 
The third possibility is to assume that UHECRs are some decay products of metastable 
superheavy particles |inpiHll2lll3lll4|ll5pi6| ll7j or topological defects ^Hl- The superheavy 
particles decay into quarks and gluons which initiate multi-hadron cascades through gluon 
bremstrahlung [13120112111221123121 • 

In the recent analysis [2] the authors have shown that for a fixed set of multiplets 
the minimal density of sources can be obtained by assuming a delta-function distribu- 
tion for h(j). We studied both this limiting luminosity, h{j) = 6{j — j-t), and a more 
realistic one with Schechter's luminosity function, which can be given as: h{j)dj = h ■ 
( J /j* ) ~ ^ ■ 2^ exp ( - j /j^ )d{j lU ) • 

The space distribution of sources can be given based on some particular survey of 
the distribution of nearby galaxies or on a correlation length rg characterizing the clus- 
tering features of sources. For simplicity the present analysis deals with a homogeneous 
distribution of sources. 

In order to determine the confidence level (CL) regions for the source densities we 
used the frequentist method [2S|- We wish to set limits on S, the source density. Using 
our Monte-Carlo based P{r, E, Ec) functions and our analytical technique we determined 
p{Ni, N2,N3, S;j^), which gives the probability of observing A''i singlet, doublet, N3 
triplet etc. events if the true value of the density is S and the central value of luminosity 
is j*. For a given set of {Ni,i = 1,2,...} the above probability distribution as a function 
of S and j^, determines the 68% and 95% confidence level regions in the S — j* plane. 
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Figure 2: The la (68%) and 2a (95%) confidence level regions for and the source density (14 
UHECR with one doublet U). 



Fig. 121 shows the confidence level regions for one of our models (with injected energy 
distribution c{E) oc E~^; and Schechter's luminosity distribution: h{j)dj oc (j/i*)~^'^^ • 
&xp{—j/j^)d{j/j:^)). The regions are deformed, thin ellipse-like objects. For this model our 
final answer for the density is 180^^gg^|^^^^'* • 10"'^ Mpc~^, where the first errors indicate the 
68%, the second ones in the parenthesis the 95% CLs, respectively. The choice of |2 -h(j)oc 
(5(j)- and, e.g. E~'^ energy distribution gives much smaller value: '^■'^'^'^L^^^^^^yo) Mpc~^, 
which is in a quite good agreement with the result of Ref . [2] . 



4. Decay of GUT scale particles 

An interesting idea discussed by refs. |15pi6pi7j is that SPs could be the source of UHECRs. 
Note, that any analysis of SP decay covers a much broader class of possible sources. Several 
non-conventional UHECR sources produce the same UHECR spectra as decaying SPs. We 
studied the scenario that the UHECRs are coming from decaying SPs and we determined 
the mass of this X particle mx by a detailed analysis of the observed UHECR spectrum 

The hadronic decay of SPs yields protons. They are characterized by the fragmentation 
function (FF) D{x,Q'^) which gives the number of produced protons with momentum 
fraction x at energy scale Q. For the proton's FF at present accelerator energies we used 
ref. [2n|- We evolved the FFs in ordinary and in super symmetric QCD to the energies of 
the SPs. This result can be combined with the prediction of the MLLA technique, which 
gives the initial spectrum of UHECRs at the energy mx (cf. Fig.. EJ. Similar results are 
obtained by P7] . 

Depending on the location of the source -halo or extragalactic (EG)- and the model 
-SM or MSSM- we studied four different scenarios. In the EG case protons loose some 
fraction of their energies, described by P{r, E, Ec). We compared the predicted and the 
observed spectrums by a maximum likelihood analysis. This analysis gives the mass of the 
SP and the error on it. 
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Figure 3: The quark FFs at Q=10^^ GeV for proton/pion in SM (solid/dotted line) and in MSSM 
(dashed/dashed-dotted line) 0. We change from logarithmic scale to linear at a; = 0.01. 
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Figure 4: UHECR data with their error bars and the best fit from a decaying SP 0. There are 
no events above 3 x lO^'^ eV (shown by an arrow). Zero event does not mean zero flux, but an 
upper bound for the flux. Thus, the experimental flux is in the "hatched" region with 68% CL. 



Fig. 0] shows the measured UHECR spectrum and the best fit, which is obtained in 
the EG-MSSM scenario. 

To determine the most probable value for the mass of the SP we studied 4 scenarios. 
Fig. [3 contains the Xmm values and the most probable masses with their errors for these 
scenarios. 

The UHECR data favors the EG-MSSM scenario. The goodnesses of the fits for the 
halo models are far worse. 

The SM and MSSM cases do not differ significantly. The most important message is 
that the masses of the best fits (EG cases) are compatible within the error bars with the 
MSSM gauge coupling unification GUT scale: mx = 10^ GeV, where b = lA.6't\'^- 
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Figure 5: The most probable values for the mass of the decaying ultra heavy dark matter with 
their error bars and the total values [Sj. Note that 21 bins contain nonzero number of events 
and the fit has 3 free parameters. 



5. Z-burst scenario 

Already in the early eighties there were some discussions about the possibility that the 
ultrahigh energy neutrino spectrum could have absorption dips at energies around E^^^ = 
M^/(2m^,.) = 4.2 • 10^^ (1 eV/m,y.) eV due to resonant annihilation with relic neutrinos 
of mass ruui ' Predicted by the hot Big Bang, into Z bosons of mass Mz [IHllSn] . Recently 
it was realized that the same annihilation mechanism gives a possible solution to the GZK 
problem It was argued that the UHECRs above the GZK cutoff are from these Z- 
bursts. The Z-burst hypothesis for the ultrahigh energy cosmic rays was discussed in many 

We compared this scenario with observations [Bj. 

The density distribution of Rz^s as hot dark matter follows the total mass distribution; 
however, it is less clustered. To take this into account, the shape of the n^^ir) distribution 
was varied, for distances below 100 Mpc, between the standard cosmological homogeneous 
case and that of the total mass distribution obtained from peculiar velocity measurements. 
These measurements suggest relative overdensities of at most a factor fi^ = 2-^3, depending 
on the grid spacing. A relative overdensity f,y = 10^ -r- in our neighbourhood, as it 
was assumed in earlier investigations of the Z-burst hypothesis, seems unlikely in view 
of these data. Note, that our quantitative results turned out to be rather insensitive to 
the variations of the overdensities within the considered range. For scales larger than 100 
Mpc the relic neutrino density was taken according to the big bang cosmology prediction, 
n^,. = 56 • (1 -|- zY cm~^. 

We gave the energy distribution of the produced particles in our lab system, which is 
obtained by Lorentz transforming the CM collider results. We included in our analysis the 
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protons, which are directly produced in the Z-burst and appear as decay products of the 
neutrons. Photons were also taken into account. They are produced in hadronic Z decays 
via fragmentation into neutral pions, Z^7r^ + X^2^ + X. Electrons (and positrons) 
from hadronic Z decay are also relevant for the development of electromagnetic cascades. 
They stem from decays of secondary charged pions, Z vr^ -\- X + X. 

The UHECR flux from Z-bursts is proportional to the differential fluxes Fi,. of ultrahigh 
energy cosmic neutrinos. Unfortunately, the value of these fluxes is essentially unknown. 
In this situation of insufficient knowledge, we took the following approach concerning the 
flux of ultrahigh energy cosmic neutrinos of type i, Fy^{Ey^,r). It was assumed to have the 
form 



where z is the redshift and where a characterizes the cosmological source evolution. Note, 
however, that, independently of the production mechanism, neutrino oscillations result in 
a uniform mixture for the different mass eigenstates. 

The next ingredient of our analysis is the propagation of the protons and photons from 
cosmological distances. This propagation was described by the appropriate P{r, Ep, E) 
probability functions (see Section 2). 

Finally, we compared the predicted and observed spectrum and extract the mass of 
the relic v and the necessary UHEz/ flux by a maximum likelihood analysis. Qualitatively, 
our analysis can be understood as follows. In the Z-burst scenario small relic i' mass needs 
large E^^^ to produce a Z. Large Ej,^^ results in a large Lorentz boost, thus large proton 
energy. In this way the detected energy determines the mass of the relic v. The analysis 
is completely analogous to that of the previous section. The observed flux is a sum of two 
terms, namely the flux from Z-bursts and a conventional part with power-law behaviour in 
the energy. This power-law part might be produced in our galaxy (halo model) or it might 
be produced extragalactically (EG model). 

The Z-burst determination of the neutrino mass seems reasonably robust. Fig. El shows 
the summary of our relic neutrino mass determination. For a wide range of cosmological 
source evolution (a = — 3 3), Hubble parameters h = 0.61 0.9, Qm, ^a, -Zmax = 2-^5, 
for variations of the possible relic neutrino overdensity in our GZK zone and for different 
assumptions about the diffuse extragalactic photon background, the results remain within 
the above error bars. The main uncertainties concerning the central values originate from 
the different assumptions about the background of ordinary cosmic rays. 

In the case that the ordinary cosmic rays above 10^^'^ eV are protons and originate 
from a region within the GZK zone of about 50 Mpc ("halo"), the required mass of the 
heaviest neutrino seems to lie between 2.1 eV< rrii, <6.7 eV at the 68% C.L. {a < 0), 
if we take into account the variations between the minimal and moderate universal radio 
background cases and the strong UIIE7 attenuation case. 

The much more plausible assumption that the ordinary cosmic rays above 10^^'^ eV are 
protons of extragalactic origin leads to a required neutrino mass of 0.08 eV< <1.3 eV at 
the 68 % C.L. (a < 0). In this case the predicted mass has a relatively strong dependence on 
the value of the universal radio background. Physically it is easy to understand the reason. 
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Figure 6: Summary of the masses of the heaviest neutrino required in the Z-burst scenario, with 
their 1 a (sohd) and 2 a (dotted) error bars, for the case of an extragalactic and a halo background of 
ordinary cosmic rays and for various assumptions about the diffuse extragalactic photon background 
in the radio band {a ^ 0, h ~ 0.71, 51 a/ = 0.3, 51a = 0.7, z„iax = 2 see Ref. 0). From left: strong 7 
attenuation, moderate and minimal universal radio background. 

The small radio background leads to a relatively large UHE7 fraction in the observations. 
They do not loose that much energy. Thus, smaller incoming UHEz^ energy and larger nii, 
is needed to describe the data. 

We performed a Monte Carlo analysis studying higher statistics. In the near future, 
the Pierre Auger Observatory will provide a ten times higher statistics, which reduces the 
error bars in the neutrino mass to about one third of their present values. 

Let us consider in detail the 7 ray spectra from Z-bursts, notably in the ~ 100 GeV 
region. As Fig. [7| illustrates, the EGRET measurements of the diffuse 7 background in the 
energy range between 30 MeV and 100 GeV |13] gives a strong constraint on the evolution 
parameter a. The high energy spectrum, and thus the neutrino mass, is independent of 
a, at low energies only a<0 seems to be compatible with the EGRET measurements, 
quite independently of different assumptions about the universal radio background (URB). 
These numerical findings are in fairly good agreement with other recent simulations |38j . 
These photon fluxes do not contain any contribution from direct photons emitted by the 
UHEz/ sources. As they are already close to the EGRET limit, one needs special sources 
that do not give contribution to the EGRET region. 

It should be stressed that, besides the neutrino mass, the UHEz^ flux at the resonance 
energy is one of the most robust predictions of the Z-burst scenario which can be verified or 
falsified in the near future. The required flux of ultrahigh energy cosmic neutrinos near the 
resonant energy should be detected in the near future by AMANDA, RICE, and the Pierre 
Auger Observatory, otherwise the Z-burst scenario will be ruled out (cf. Fig. If such 
tremendous fluxes of ultrahigh energy neutrinos are indeed found, one has to deal with 
the challenge to explain their origin. It is fair to say, that at the moment no convincing 
astrophysical sources are known which meet the requirements for the Z-burst hypothesis, 
i.e. which have no or a negative cosmological evolution, accelerate protons at least up to 
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Figure 7: The available UHECR data with their error bars and the best fit from Z-bursts, for 
various cosmological evolution parameters a and an energy attenuation of photons exploiting a 
"minimal" URB (h = 0.71, f^M = 0.3, f^A = 0.7, z^ax = 2 see Ref. Also shown is the diffuse 7 
background in the energy range between 30 MeV and 100 GeV as measured by EGRET (solid). 



10 eV, are opaque to primary nucleons and emit secondary photons only in the sub-MeV 
region. It is an interesting question whether such challenging conditions can be realized in 
BL Lac objects, a class of active galactic nuclei for which some evidence of zero or negative 
cosmological evolution has been found (see Ref. |43j and references therein) and which were 
recently discussed as possible sources of the highest energy cosmic rays . 

Z.F thanks the organisers of the 26th Johns Hopkins Workshop on Current Problems 
in Particle Physics for the stimulating atmosphere. This work was partially supported 
by Hung. Sci. grants No. OTKA-T37615T34980/T29803/M37071M28413/OM-MU-708/- 
OMFB-1548. 
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